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Matrix product ansatz for Fermi fields in one dimension 
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We present an implementation of a continuous matrix product state for two-component fermions in 
one dimension. We propose a construction of variational matrices with an efficient parametrization 
that respects the translational symmetry of the problem (without being overly constraining) and 
readily meets the regularity conditions that arise from removing the ultraviolet divergences in the 
kinetic energy. We test the validity of our approach on an interacting spin-1/2 system and observe 
that the ansatz correctly predicts the ground-state magnetic properties for the attractive spin-1/2 
Fermi gas, including the phase-oscillating pair correlation function in the partially polarized regime. 


An exact description of an arbitrary many-body quan¬ 
tum state can be a formidable task, as the number of 
distinct states grows exponentially on the number of par¬ 
ticles or system volume. Nevertheless, accurate approx¬ 
imations of low-energy states are often feasible, as the 
vast majority of the quantum states in the entire Hilbert 
space need not be considered [1]. Among numerical treat¬ 
ments of one-dimensional (ID) strongly correlated prob¬ 
lems, the density matrix renormalization group (DMRG) 
algorithm [2] for lattice systems has been particularly no¬ 
table for its optimal efficiency. Moreover, the DMRG 
procedure is essentially equivalent to an optimization of 
a matrix product state (MPS) [3] — a natural and ef¬ 
fective form of parametrization of variational states for 
systems with limited entanglement due to short-range in¬ 
teractions. 

Recently, a coherent-state extension of the MPS to 
quantum fields in a ID continuum (called cMPS) [4-14] 
has been pioneered, and it was successful in correctly pre¬ 
dicting the ground-state energies [4] and dispersion rela¬ 
tions [10] of interacting bosons. For fermionic systems, 
on the other hand, the cMPS has been tested on a ID 
relativistic model related to quantum chromodynamics 
[7], but its success for nonrelativistic fermions — which 
have a wide range of applications in condensed matter, 
quantum information, as well as cold-atom systems - 
has not been demonstrated thus far. In particular, as re¬ 
cent experimental techniques [15] have enabled the phys¬ 
ical realization of ultracold ID Fermi gas systems [16-19] 
with a tunable interaction, particle density, and spin po¬ 
larization, a cMPS ansatz that can reliably investigate 
those systems would be of significant interest. 

This Rapid Communication extends the previous de¬ 
velopments and proposes an implementation of a cMPS 
ansatz for two-component fermions. Although some im¬ 
portant insights on the structure of a fermionic ansatz 
have been noted previously [7-9] , a subsequent work that 
fully develops a more efficient and explicit scheme that 
i) readily satisfies the regularity conditions for a finite 
kinetic energy, as well as ii) being able to capture spon¬ 
taneous translational symmetry breaking properties, is 
much needed to advance the subject further. We shall 


show that our cMPS ansatz is able to faithfully portray 
the ground-state magnetic properties of a spin-1/2 Fermi 
gas described by the Gaudin-Yang Hamiltonian [20], in 
agreement with the exact solutions obtained from the 
thermodynamic Bethe ansatz [21, 22], including the dis¬ 
tinct phases for the attractive system and spontaneous 
translational symmetry breaking properties such as inho¬ 
mogeneous (oscillatory-phase) superconductivity in the 
partially polarized regime. 

A cMPS for spin-1/2 fermions (see Fig. 1), for a system 
length of L and periodic boundary conditions, has the 
general form [7] 

|x) = Tr aux ['Pe^° i da: [2 (:c) ® 1 +Y„ G)]j |q) ( (!) 

where Q(x), TZ a (x) G C DxD and act on a D-dimensional 
auxiliary space, D is called the bond dimension, a is the 
index for spins f and J,, / is the identity operator on 
the Fock space, Tr aux is a trace over the auxiliary space, 
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FIG. 1. (Color online) Artistic rendition of the cMPS in 
Eq. (1), that when explicitly expanded, can alternatively be 
understood as a superposition of position eigenstates, each 
state having some fixed number of particles that ranges from 
1 to oo. The probability amplitude for the particular po¬ 
sition eigenstate illustrated above for xi -1 = Xi is given 
by ■•■) = Tr[- ■ ■ lZf(xi- 2 ) 

u(xi- 2 ,Xi-i)lZ^(xi-i)Ui(xi)u(xi,x i+ i)Ui(x i+ i) ■■■], where 
u(x,y) = Vex p [JJ dxQ{x) 1 \ can be viewed as a free propa¬ 
gator and IZa(x) as a vertex that creates a particle of type o 
at position x [9]. 
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■Pexp is a path-ordered exponential, and |ff) is the Fock- 
space vacuum state. The Fock-space creation and annihi¬ 
lation field operators obey the anticommutation relation 
{ip a {x),i)l,{x')} = 6 aa ’5(x - x'). 

In order to illustrate our implementation of the cMPS 
ansatz, we shall consider the Gaudin-Yang Hamiltonian 
[20], where spin-1/2 fermions interact via a contact po¬ 
tential and can be written as 

Hgy = [ dx y 2 (dx^idx^a - 2A$^t^Vv) ,( 2 ) 
J ° <r=U 


Rf(x) + TZi(x) 0 R\,(x) = T remains translationally in¬ 
variant (the bars denote complex conjugation of matrix 
entries). Some of the correlators and observables that we 
will use below are 

(i>t'(x + Sx)ip l j(x)) = Tr [e T( - L ~ Sx \R a 0 I)e fSx (I 0 R a ')\ 

xe i(Qa-g a i)x e -iq a iSx 

{d x rp] T {x)d x ^ <T {x)) = Ti {e TL [(iq a R a + [Q, R a }) 0 c.c.]}, 
C pail (6x) = T 'r[e T ( L - Sx \R t R i 0 I)e T5x 

x(I®R t R l )}e- i ^+u'> Sx , (5) 


where a denotes the conjugate spin to a, and the in¬ 
teraction strength is A < 0 and A > 0 for repulsive 
and attractive interactions, respectively. The exact wave 
function for the two-particle spin-singlet sector of this 
translationally invariant system with periodic boundary 
conditions can be obtained using the coordinate Bethe 
ansatz as [22] 

(f>s (. XX , x 2 ) = A s e^ kl+k2 ^ Xl+X2) 

x + Sse -i(k 2 -kt)\x\\ ^ 

where kj is the quasimomentum of the jth particle (obey¬ 
ing the corresponding Bethe ansatz equation), x = x 2 — 
xi, S s = (fa — k 2 — 2i A) 1 (fa — k 2 + 2iA), and A s is 
an overall amplitude. By constructing a cMPS with 


Q{x) 

U t {x) 

Kl(x) 
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-diag (q- -q + k 2 - fa fa - k 2 
e iq + x f a+ 2 1 2 \ 

2 ^ 0 -(7+ ) ’ 

A s (l + S)e iq ~ x fa- 0 \ 

4isin[(£: 2 — fa)L/2\ \ —2I 2 —a _ J 


q+ - q-), 
( 4 ) 


where fa is a d x d identity matrix, a± = a x ± ia y with 
Pauli matrices cr x , y ,z> and q± = (fa + k 2 ) /2±mr/L with 
n £ Z (its parity a function of k 2 — fa), we recover the 
exact wave function in Eq. (3) after projecting the cMPS 
onto the two-particle sector. 

We remark the phase modulation present on the lZ a 
matrices in Eq. (4), while the Hamiltonian is transla¬ 
tionally invariant, is reminiscent of the general result for 
Bloch (or Floquet) states. In fact, it can be shown that 
e i q„x j g mos t general form of modulation on other¬ 
wise (spatially) constant TZ a matrices that preserves the 
continuous translational symmetry of the system. Based 
on this insight, we depart from the extended practice 
of taking both Q(x) and lZ a (x) to be independent of 
coordinates, the so-called uniform ansatz or uMPS [23], 
and propose a phase modulation on 1Z„(x). Henceforth, 
Q(x) = Q and R a (x) = R a e lq<,x , where Q and R a are 
D x D matrices that are independent of position, and 
q a are real variational parameters, unconstrained in the 
thermodynamic limit. We note that despite this modu¬ 
lation of Tier, T(x) = Q(x) 0 I + I 0 Q(x) + 7&|-(x) 0 


where T = Q®I+I®Q — i?-|-0i?f — R±®R±, C pa i r (Sx) = 
(ip\(x + 5x)ip\(x + Sx)ipi(x)ip^(x)), and 5x > 0. These 
expressions are independent of x (with a' = a in the 
first case) and, in general, translationally invariant ex¬ 
pectation values are expected from our formalism when¬ 
ever the operator conserves the total particle number and 
spin. Likewise, similar expressions and conclusions can 
be drawn for Sx < 0. 

The kinetic energy density, written in terms of q a , Q , 
and R a , can be derived from a lattice-regularized expec¬ 
tation value ({e _1 [?/4 (x n+ 1 ) — (x„)]} • {h.c.}) where e 

= x n +i — x n . It contains terms having powers of e -1 and 
e~ 2 that diverge in the continuum limit of e —> 0. These 
divergent terms cancel out exactly in the bosonic cMPS, 
whereas the kinetic energy density expression needs to 
be regularized [9] in the fermionic cMPS. The regular ex¬ 
pression for the kinetic energy density and the simplified 
form for C pa i r (<5x), both stated in Eqs. (5), follow after 
imposing the conditions: 

{R'l, Rfa = 0 and R 2 = 0. (6) 

The optimal method to meet these conditions is to con¬ 
struct Rf and R± to have special algebraic structures. 
One systematic way to enforce the requirement R 2 = 0 
is by having R^ in a Jordan canonical form [24] with 
only nontrivial Jordan blocks. As the nilpotency degree 
of R -j- must be 2, all of its Jordan blocks must be 2 x 2, 
and we are lead to the unique choice Rf = Id / 2 0 er+/ 2 - 
A general form for R^ that satisfies both R 2 = 0 and 
{R t ,Rfa = 0 is 

R i = [- P d/2( 7 o/ 4 ® 0-+/2)Po/ 2 ] ® cr~ (7) 

+ [-^0/2 (^W 4 ® CT +/ 2 + R D/4 ® fa)PD/2 ] ® 

where A D/4 , B D/4 <E C D / 4xD / 4 , and P D/2 € C D / 2xD / 2 
is an invertible matrix. 

Furthermore, we enforce Q + + R^R^ + = 0 

using the cMPS gauge freedom [4, 9] , so that the general 
form of Q becomes Q = A — — ^R^R^, where A is 

an anti-Hermitian matrix that is independent of position. 
Interestingly, in our construction of the Q and R a ma¬ 
trices, the R^ matrix contains no variational parameters, 
and all parameters reside on the Q and matrices. All 
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FIG. 2. (Color online) Zero temperature ground-state phase 
diagram for the homogeneous attractive spin-1/2 Fermi gas 
in ID. The phase boundaries were obtained from the Bethe 
ansatz [26]. The points X, Y, and Q along the circular tra¬ 
jectory C denote the intersections with the phase boundaries. 


FIG. 3. (Color online) Free energy densities along the trajec¬ 
tory C in Fig. 2 obtained from variational states having D = 
4,8,16,20 and the exact Bethe ansatz result (solid curve). 
The inset has free energies from variational states (D = 16) 
with lZ a (x) = R a (solid circles) and TZ a {x) = R a e tq,7X (open 
circles). 


the information needed to correctly capture the physics 
of the spin-up particles is, nevertheless, correctly encoded 
in Q and the transfer matrix T.[25] 

We note that our form of R a naturally constrains the 
bond dimension D of the cMPS to be a multiple of 4, 
and the total number of variational degrees of freedom 
between the Q and R a matrices becomes approximately 
1.75D 2 . For a given bond dimension D 1 it is desirable to 
retain the maximum number of nonredundant variational 
parameters, as the computational time required to com¬ 
pute the expectation values for an arbitrary state scales 
as 0(D 6 ). Our implementation gives an improvement 
of nearly a factor of 2 as compared to a previous con¬ 
struction that had D 2 + D variational degrees of freedom 
[ 8 ]. 

In order to test our ansatz, we refer again to the 
Gaudin-Yang system in a grand canonical ensemble un¬ 
der the presence of a constant Zeeman field that is 
perpendicular to the system axis. Figure 2 shows the 
ground-state phase diagram for the attractive Fermi gas 
in ID, which has four distinct phases that are character¬ 
ized by the population of each of the fermions: (i) fully 
polarized when n-j- > = 0 (between X and Y along the 

trajectory C), (ii) partially polarized when n-j- > nj, > 0 
(between Y and *P), (iii) fully paired when = nj. > 0 
(between d' and D), and (iv) vacuum when n-f- = nj. = 0 
(between f l and X). 

To obtain variational ground states at the various 
phases, we optimized the expectation value of the zero 
temperature free energy density / = F/L in the thermo¬ 
dynamic limit. For the minimization, we implemented 
the simulated annealing algorithm [27], taken as a variant 
of the simplex method that incorporates a Metropolis¬ 
like scheme [28] . The free energy density operator is / = 
H G y/L-h (n t + %)-/i (n t - %), where n = (/x t +/x.0/2 


is the chemical potential, h = (/if — /q,)/2 the effective 
magnetic field, and n a the densities of the two species. In 
terms of the variational matrices, we evaluated (/) using 
the expressions in (5). 

The free energy densities obtained from the optimiza¬ 
tion of the cMPS variational ansatz along the trajectory 
C are plotted in Fig. 3. The bond dimensions range from 
D = 4 to 20. While a qualitative agreement starts to 
emerge from D > 8, the variational ansatz with D > 16 
yields quantitatively accurate approximations, and con¬ 
vergence to the exact results is readily seen as D is fur¬ 
ther increased. The inset in Fig. 3 presents the effect of 
neglecting the spatial modulation on the lZ a {x) matri¬ 
ces by enforcing <7f = q± = 0. In the partially polarized 
regime - i.e. between the points Y and 'P — we see 
notable gaps in free energy densities between the two 
curves, an indication that the modulation of lZ a (x) is 
crucial for the accurate approximation of ground states 
in the partially polarized phase. Moreover, we find that 
the relation q^/qi = —ni/rif holds true in the partially 
polarized regime. 

Figure 4 is a plot of the densities nf, n f and their 
difference along the trajectory C. As in the free energy 
density plot in Fig. 3, D = 8 only gives results that qual¬ 
itatively agree with the exact results, while D > 16 show 
much more quantitatively reliable approximations. The 
predictions for the phase boundary points X, Y and D 
are excellent, whereas the transition point *P between the 
partially polarized and the fully paired phases converges 
more slowly; gradual improvement of the prediction with 
increasing D is evident, nonetheless. 

One of the exciting recent investigations by the 
ultracold-atom community is the effort to experimentally 
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FIG. 4. (Color online) Densities of the majority spins n-[ 
(blue), minority spins rq, (red), and their difference n-f — nj. 
(green) along the trajectory C in Fig. 2, as obtained from 
variational states having bond dimensions D = 8,16, 20 and 
the exact Bethe ansatz result (solid curve). 


confirm whether the ground state of a spin-imbalanced 
spin-1/2 Fermi superfluid is, in fact, the Fulde-Ferrell- 
Larkin-Ovchinnikov (FFLO) state [29]. This state, pos¬ 
tulated about 50 years ago, is characterized by a pair 
correlation function ('ipUx + 6x)^r,(x + 
that has a long-range oscillatory behavior. The renewed 
excitement was stirred by a theoretical prediction that a 
large portion of the ground-state phase diagram for an 
attractive ID spin-1/2 Fermi gas is the FFLO phase (as 
seen in Fig. 2) [26, 30-33], paving an avenue for a feasible 
experimental verification [16, 19]. 

One of the advantages of an MPS formalism is its abil¬ 
ity to express expectation values of essential correlations 
and observables in concise, closed forms. Using the cMPS 
approximation of the ground state at point A in Fig. 2, we 
have evaluated the pair correlation function in the par¬ 
tially polarized regime via Eq. (5). In Fig. 5, an oscilla¬ 
tion that is present over a long distance is clearly visible, 
agreeing with the prediction that the ground state real¬ 
izes the FFLO physics. We note that this behavior was 
predicted beforehand, as e TSx converges to a fixed value 
for large Sx and the only remaining spatial dependence 
becomes On the other hand, the fully paired 

phase does not show an oscillatory pair correlation, and 
we verified there is no long-range correlation for the same 
system with a repulsive interaction. 

We have shown the validity of our implementation of 
a fermionic cMPS that incorporates a spatial modula¬ 
tion on the lZ a matrices with an efficient, explicit alge¬ 
braic construction to meet the regularity conditions for 
the kinetic energy. We have verified that our fermionic 
cMPS yields accurate results for interacting spin-1/2 
Fermi gases in ID (that converge to the exact results ob¬ 
tained from the thermodynamic Bethe ansatz), and have 



FIG. 5. (Color online) Real part of the pair correlation 
function in the partially polarized phase (black) and the fully 
paired phase (red), evaluated from variational cMPS (D = 16) 
at points A and B respectively (see Fig. 2). The black dashed 
curve is the imaginary part of the pair correlation at A. The 
blue curve is obtained from a D = 16 cMPS in a repulsive 
system, with A = —1, fj,/ |A| 2 = 0.5, and h/\A\ 2 = 0.01. All 
correlators have been normalized by n 2 , where n = n-f + n_\,. 
Inset: Fourier transform n a {k) of ('/>J.(:r)'i/v(0)) evaluated at 
points A and B, normalized by n. At large k , n-ik) decays 
as k~ 4 (cf. Refs. [7, 34]). 

found that the variational states correctly predict essen¬ 
tial observables and correlations, including a spontaneous 
translational symmetry breaking phenomenon, such as 
the FFLO superconductivity. As the cMPS ansatz is 
generic and not restricted to integrable Hamiltonians, we 
envisage that our work would serve as a critical stepping 
stone towards solving outstanding ID problems and shed 
light on our understanding of interacting Fermi fields. 
Future study includes the extension of the cMPS to vari¬ 
ous interesting condensed-matter and cold-atom systems, 
such as spin-orbit coupled nanowires [35], FFLO super¬ 
fluids in traps [33, 34, 36] or tubular lattices [30, 37], and 
Bose-Fermi mixtures [38]. 
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ARO W911NF-07-1-0464; parallel computing resources 
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